The present invention relates to a method for adaptively estimating, with a projection algorithm, a transfer function of an unknown system and its output in an acoustic canceller, active noise control or the like and an estimating device using such a method.
In the following description, time will be represented by a discrete time k. For example, the amplitude of a signal x at time k will be expressed by x(k). FIG. 1 is a diagram for explaining the estimation of the transfer function of an unknown system. Reference numeral 11 denotes a transfer function estimation part and 12 the unknown system, and reference character x(k) represents an input signal to the unknown system and y(k) an output signal therefrom. The transfer function h(k) of the unknown system is estimated using the input signal x(k) and the output signal y(k). FIG. 2 is a diagram for explaining an adaptive estimation of the transfer function. Reference numeral 21 denotes an estimated transfer function correcting vector calculation part, 22 an estimated transfer function correction part and 23 a convolution part, these parts constituting an estimated signal generation part 20. The transfer function h(k) of the unknown system 12 is estimated as a transfer function h(k) of an FIR filter of a tap number L which forms the convolution part 23. More specifically, coefficients h.sub.1 (k) . . . h.sub.L (k) of the FIR filter are estimated. Let it be assumed that the "transfer function" and the "FIR filter coefficient" will hereinafter be construed as the same. For the sake of brevity, the filter coefficient is represented as an estimated transfer function vector h(k) defined by the following equation. EQU h(k)=[h.sub.1 (k), h.sub.2 (k), . . . , h.sub.L (k)].sup.T ( 1)
where T represents a transpose.
In FIG. 2, the input signal x(k) to the unknown system 12 is fed to the convolution part 23 and a calculation is performed to obtain the estimated transfer function vector h(k) that minimizes an expected value which is the square of an error signal e(k) available from a subtractor 24 which detects a difference between an output y(k) from the convolution part 23 given by the following equation (2) and the output y(k) from the unknown system 12. EQU y(k)=h(k).sup.T x(k) (2) EQU x(k)=[x(k), x(k-1), . . . , x(k-L+1)].sup.T ( 3)
where y(k) is an estimated value of the output from the unknown system, which is close to the value of the output y(k) when the estimated transfer function h(k) is close to an unknown characteristic.
In practice, the transfer function of an unknown system often varies with time as in the case where the transfer function of an acoustic path varies with movement of audiences or objects in a sound field such as a conference hall or theater. For this reason, the estimated transfer function h(k) of the unknown system also needs to be adaptively corrected in accordance with a change in the acoustic path of the unknown system. The estimated transfer function correcting vector calculation part 21 calculates an estimated transfer function correcting vector .delta.h(k) on the basis of the error signal e(k) and the input signal x(k) to the unknown system 12. The estimated transfer function correction part 22 corrects the estimated transfer function by adding the estimated transfer function correcting vector .delta.h(k) to the estimated transfer function vector h(k) as expressed by the following equation. EQU h(k+1)=h(k)+.mu..delta.h(k) (4)
where .mu. is called a step size, which is a preselected quantity for controlling the range of each correction and is handled as a time-invariant constant. In the following description, the estimated transfer function correcting vector .delta.h(k) is calculated for .mu.=1 and, if necessary, it is multiplied by a desired step size .mu.. In some applications the value of the step size .mu. is caused to vary with time, but in such a case, too, the following description is applicable. The corrected estimated transfer function vector is transferred to the convolution part 23. The above is a transfer function estimating operation at time k and the same operation is repeated after time k+1 as well.
The estimated transfer function correcting method described above in respect of FIG. 2 is known as an adaptive algorithm. While an LMS (Least Mean Square) algorithm and an NLMS (Normalized LMS) are also well-known as adaptive algorithms, a description will be given of a projection algorithm proposed in a literature [Ozeki and Umeda: "An Adaptive Filtering Algorithm Using an Orthogonal Projection to an Affine Subspace and Its Properties", Journal of Institute of Electronics, Information and Communication Engineers of Japan (A), J67-App. 126-132, (1984-2).]
The projection algorithm requires a larger number of operations than does the NLMS but has an excellent adaptability to a speech signal input. With the projection algorithm, as referred to above, the vector .delta.h(k) is determined by Eq. (4) for .mu.=1 so that simultaneous equations composed of the following p equations are satisfied. ##EQU2## Eq. (5) indicates that the vector .delta.h(k) is determined so that the estimated transfer function h(k+1) updated at time k provides the same values y(k), y(k-1), . . . , y(k-p+1) as the outputs from the unknown system, respectively, for p input vectors x(k), x(k-1), . . . , x(k-p+1) prior to time k. By this, it is expected that the characteristic of the estimated transfer function h(k) will approach the characteristic of the unknown system as the adaptive updating of the estimated transfer function is repeated using the vector .delta.h(k).
In the above, p is a quantity commonly called a projection order. As the projection order p increases, the adaptability of the projection algorithm increases but the computational complexity also increases. The conventional NLMS method corresponds to the case of p=1.
Now, transposing the first equation in Eq. (5), we have EQU x(k).sup.T .delta.h(k)=y(k)-x(k).sup.T h(k)=e(k) (6)
Furthermore, transposing the second equation in Eq. (5) and using an equation obtainable by setting k in Eqs. (4) and (6) to k-1, we have ##EQU3## Thereafter, the following relationship similarly holds. EQU x(k-i).sup.T .delta.h(k)=(1-.mu.).sup.i e(k-i) (8)
Based on this, Eq. (5) can be rewritten as the following system of simultaneous equations. EQU X.sub.p (k).delta.h(k)=e(k) (9)
where X.sub.p (k) is a matrix with p rows and L columns and e(k) is a vector of the p rows; they are defined by the following equations. ##EQU4## The vector e(k) will hereinafter be referred to as an error vector. Now, since p is usually smaller than L, Eq. (9) is an under-determined simultaneous equation or indeterminate equation for the vector .delta.h(k) and the minimum norm solution of the vector .delta.h(k) is given by the following equation. ##EQU5## where EQU g(k)=R.sub.p (k).sup.-1 e(k) (13) EQU R.sub.p (k)=X.sub.p (k)X.sub.p (k).sup.T ( 14)
R.sub.p is a matrix with p rows and p columns, which will hereinafter be referred to as a p-order covariance matrix or auto-correlation matrix, and g(k) a pre-filter coefficient vector. Letting elements of the pre-filter coefficient vector g(k) be represented by g.sub.1 (k), g.sub.2 (k), . . . , g.sub.p (k), the estimated transfer function correcting vector .delta.h(k) can be expressed on the basis of Eq. (12) as follows: ##EQU6## When, when the projection algorithm is used in accordance with the present invention, the estimated transfer function correcting vector calculation part 21 in FIG. 2 has such a construction as depicted in FIG. 3. Reference numeral 31 denotes a pre-filter coefficient vector calculation part, which uses the input signal x(k) and the error signal e(k) to calculate the pre-filter coefficient g(k) on the basis of Eq. (13). Reference numeral 32 denotes a pre-filtering part, which performs the pre-filtering operation expressed by Eq. (15) to synthesize the estimated transfer correcting vector .delta.h(k) by use of the pre-filter coefficient g(k) that is transferred from the pre-filter coefficient vector calculation part 31.
Now, a description will be given of the computational complexity of the projection scheme described above. The computational complexity mentioned herein is the number of multiplication-addition (or addition) operations necessary for estimating operations per unit discrete time. The computational complexity of Eq. (2) in the convolution part 23 of the tap number L in FIG. 2 is L. The computational complexity of Eq. (13) in the pre-filter coefficient vector calculation part 31 of the estimated transfer function correcting vector calculation part 21 is about p.sup.3 /6 when using the Choleski method which is a typical computation method. The computational complexity of Eq. (15) in the pre-filtering part 32 is (p-1)L. The computational complexity of Eq. (4) in the estimated transfer function correction part 22 is L. Thus, the entire computational complexity NC of the projection scheme is given by the following equation. EQU NC=L+p.sup.3 /6+(p-1)L+L (16)
On the other hand, the computational complexity of the NLMS scheme or LMS algorithm is about 2L. For example, when L=500 and p=20 (a typical value in the case of an acoustic echo canceler), the number of operations involved in the NLMS scheme is 1000, whereas the projection scheme requires as many as about 12000 operations on the basis of Eq. (16). The computational complexity p.sup.3 /6 of Eq. (13), in particular, abruptly increases as the projection order p increases. Thus, the projection scheme has excellent convergence characteristics as compared with the NLMS scheme but poses the problem of increased computational complexity.